# -*- coding: utf-8 -*-
from scipy import *
from scipy import weave
from scipy.linalg import inv,eig
from scipy.sparse.linalg import eigs

Import_Mayavi = '''
try: from enthought.mayavi.mlab import *
except: from mayavi.mlab import *
'''
Import_Pylab='''
from matplotlib.pyplot import *
rcParams['lines.linewidth']=3.
rcParams['lines.color']='k'
rcParams['lines.markersize']=8
rcParams['lines.markeredgewidth']=1.5
rcParams['legend.fontsize']=25
rcParams['font.size']=28
rcParams['font.family']='serif'
rcParams['font.serif']='Times New Roman'
rcParams['axes.formatter.limits']=[-3,4]
'''
D = lambda E,t,v: E*t**3/12./(1-v**2)
delta = lambda p,q: 1. if p==q else 0.
B1 = lambda p,q: 0.5 if p==q else 0.
B2 = lambda p,q: 0. if (p+q)%2==0 else p*2./(p**2-q**2)

H1 = lambda i,j,m,n: (i+j)**2*(m+n)**2*B1(i+j,m+n)
H2 = lambda i,j,m,n,beta: (i*i-beta*beta*j*j+i*j*(1-beta*beta))*(m*m-beta*beta*n*n+m*n*(1-beta*beta))*B1(i+j,m+n)

def solve_eigen(U,T,n_modes):

	K,A=eigs(dot(inv(U),T),n_modes)
	K,A = 1./K.real,A.real
	sgn = {}
	KA = {}
	for i in range(len(K)):
		if abs(K[i]) in KA.keys():
			KA_key = abs(K[i])+1.0e-10
		else:
			KA_key = abs(K[i])
		KA[KA_key] = A[:,i]
		sgn[KA_key] = K[i]/abs(K[i])
		
	K = array(sorted(KA.keys()))
	A = array(map(KA.get, K)).T
	K = K*array(map(sgn.get, K))
	
	return K,A
	
def platebuckling_sim(M,N,beta=1.0,gamma=0.0,n_modes=6,ds=0):
	'''
	Elastic Buckling Analysis of four-edge simple supported rectangular plate with diagnal stiffeners.
	Input: 
		M,N: Number of series to approximate transverse deflection of the plate;
		beta: aspect ratio of the plate;
		gamma: stiffenning ratio;
		n_modes: Number of buckling modes to calculate;
	Output:
		K: eigen values;
		A: factors of series
	'''
	
	Nut = M*N
	U,T = zeros((Nut,Nut)),zeros((Nut,Nut))
	
	for P in range(Nut):
		for Q in range(P,Nut):
			i,j = P/N+1,P%N+1
			m,n = Q/N+1,Q%N+1
			U[P,Q] = beta*pi**2/4.*((i/beta)**2+j**2)**2*delta(m,i)*delta(n,j)+0.25*gamma*(beta**2+1.)**(-1.5)*pi**2*(1.0+ds*(-1.)**(j+n))*(H1(i,j,m,n)-H1(i,j,m,-n)-H1(i,-j,m,n)+H1(i,-j,m,-n))
			T[P,Q] = -2.*i*j*B2(m,i)*B2(n,j)
	
	U = triu(U,1).T+U
	T = triu(T,1).T+T
	
	K,A=solve_eigen(U,T,n_modes)

	return K,A

def platebuckling_fix(M,N,beta=1.0,gamma=0.0,n_modes=6,ds=0):
	'''
	Elastic Buckling Analysis of four-edge clapped rectangular plate with diagnal stiffeners.
	Input: 
		M,N: Number of series to approximate transverse deflection of the plate;
		beta: aspect ratio of the plate;
		gamma: stiffenning ratio;
		n_modes: Number of buckling modes to calculate;
	Output:
		K: eigen values;
		A: factors of series
	'''
	
	Nut = M*N
	U,T = zeros((Nut,Nut)),zeros((Nut,Nut))
	
	for P in range(Nut):
		for Q in range(P,Nut):
			i,j = P/N,P%N
			m,n = Q/N,Q%N
			U[P,Q] = beta*pi**2/4.*((i/beta)**2+j**2)**2*delta(m,i)*delta(n,j)*(1+delta(i,0)+delta(0,j))+0.25*gamma*(beta**2+1.)**(-1.5)*pi**2*(1.+ds*(-1.)**(j+n))*(H1(i,j,m,n)+H1(i,j,m,-n)+H1(i,-j,m,n)+H1(i,-j,m,-n))
			T[P,Q] = 0. if i*j*m*n==0 else (i**2*n*1.0/m+j**2*m*1.0/n)*B2(m,i)*B2(n,j)
	
	U = triu(U,1).T+U
	T = triu(T,1).T+T
	
	Nbc = 2*M+2*N
	BC = zeros((Nut,Nbc))
	
	for i in range(M):
		for j in range(N):
			BC[i*N+j,j] = 1.
			BC[i*N+j,N+j] = (-1.)**i
			BC[i*N+j,2*N+i] = 1.
			BC[i*N+j,2*N+M+i] = (-1.)**j

	
	K_U = vstack((hstack((U,BC)),hstack((BC.T,zeros((Nbc,Nbc))))))
	K_T = zeros((Nut+Nbc,Nut+Nbc))
	K_T[:Nut,:Nut] = T
	
	K,A=solve_eigen(K_U,K_T,n_modes)

	return K,A

def platebuckling_simfix(M,N,beta,gamma,n_modes=6,ds=0):
	'''
	Elastic Buckling Analysis of two-edge clapped rectangular plate with diagnal stiffeners.
	Input: 
		M,N: Number of series to approximate transverse deflection of the plate;
		beta: aspect ratio of the plate;
		gamma: stiffenning ratio;
		n_modes: Number of buckling modes to calculate;
	Output:
		K: eigen values;
		A: factors of series
	'''

	M,N = int(M),int(N)
	beta,gamma = float(beta),float(gamma)
	Nut = int(M*N)
	U,T = zeros((Nut,Nut)),zeros((Nut,Nut))

	Nbc = M+N
	BC = zeros((Nut,Nbc))
	
	for P in range(Nut):
		for Q in range(P,Nut):
			i,j = P/N,P%N
			m,n = Q/N,Q%N
			U[P,Q] = beta*pi**2/4.*((i/beta)**2+j**2)**2*delta(m,i)*delta(n,j)*(1.+delta(0,j))+0.25*gamma*(beta**2+1.)**(-1.5)*pi**2*(1.+ds*(-1.)**(j+n))*(H1(i,j,m,n)+H1(i,j,m,-n)+H1(i,-j,m,n)+H1(i,-j,m,-n))
			T[P,Q] = 0. if i*j*m*n==0 else (1.0*i*n+1.0*i*j*j/n)*B2(m,i)*B2(n,j)
	
	for i in range(M):
		for j in range(N):
			BC[i*N+j,i] = 1.
			BC[i*N+j,M+i] = (-1.)**j

	U = triu(U,1).T+U
	T = triu(T,1).T+T

	K_U = vstack((hstack((U,BC)),hstack((BC.T,zeros((Nbc,Nbc))))))
	K_T = zeros((Nut+Nbc,Nut+Nbc))
	K_T[:Nut,:Nut] = T

	K,A=solve_eigen(K_U,K_T,n_modes)

	return K,A
	
def platebuckling_simw(M,N,beta,gamma,n_modes=6,ds=0):
	'''
	Elastic Buckling Analysis of four-edge simple supported rectangular plate with diagnal stiffeners.
	M,N: Number of series to approximate transverse deflection of the plate;
	beta: aspect ratio of the plate;
	gamma: stiffenning ratio;
	phi: ratio;
	n_modes: Number of buckling modes to calculate;
	'''

	M,N = int(M),int(N)
	beta,gamma,ds = float(beta),float(gamma),float(ds)
	Nut = int(M*N)
	U,T = zeros((Nut,Nut)),zeros((Nut,Nut))
	
	code = '''
	int p,q,i,j,m,n;
	const float PI = 2.*acos(0.); 
	for (p = 0; p < Nut; p++)
	{   
	for (q = p; q < Nut; q++)
	{
	i = p/N+1; j = p%N+1;
	m = q/N+1; n = q%N+1;
	U(p,q) = beta*PI*PI/4.*pow((pow(i/beta,2)+j*j),2)*delta(m,i)*delta(n,j)+0.25*gamma*pow((1.+beta*beta),-1.5)*(1.+ds*pow(-1.,(n+j)))*PI*PI*(pow((i+j),2)*pow((m+n),2)*B1(i+j,m+n)-pow((i+j),2)*pow((m-n),2)*B1(i+j,m-n)-pow((i-j),2)*pow((m+n),2)*B1(i-j,m+n)+pow((i-j),2)*pow((m-n),2)*B1(i-j,m-n));
	T(p,q) = -2.0*i*j*B2(m,i)*B2(n,j);
	}
	}

	'''
	extra_code = '''
	float delta(int i, int j)
	{
	if (i == j) {return 1.;}
	else {return 0.;}
	}

	float B1(int i, int j)
	{
	if (i == j) {return 0.5;}
	else {return 0.;}
	}

	float B2(int i, int j)
	{
	if ((i+j)%2==0) {return 0.;}
	else {return i*2./(pow(i,2)-pow(j,2));}
	}

	'''
	weave.inline(code, ['U','T','Nut','N','beta','gamma','ds'],
	type_converters=weave.converters.blitz,
	support_code=extra_code,
	headers = ['<math.h>'],
	compiler = 'gcc')

	U = triu(U,1).T+U
	T = triu(T,1).T+T

	K,A=solve_eigen(U,T,n_modes)

	return K,A

def platebuckling_simfixw(M,N,beta,gamma,n_modes=6,ds=0):
	'''
	Elastic Buckling Analysis of two-edge clapped rectangular plate with diagnal stiffeners.
	Input: 
		M,N: Number of series to approximate transverse deflection of the plate;
		beta: aspect ratio of the plate;
		gamma: stiffenning ratio;
		n_modes: Number of buckling modes to calculate;
	Output:
		K: eigen values;
		A: factors of series
	'''

	M,N = int(M),int(N)
	beta,gamma,ds = float(beta),float(gamma),float(ds)
	Nut = int(M*N)
	U,T = zeros((Nut,Nut)),zeros((Nut,Nut))

	Nbc = M+N
	BC = zeros((Nut,Nbc))

	code = '''
	int p,q,i,j,m,n;
	const float PI = 2.*acos(0.); 
	for (p = 0; p < Nut; p++)
	{   
	for (q = p; q < Nut; q++)
	{
		i = p/N+1; j = p%N;
		m = q/N+1; n = q%N;
		
		U(p,q) = beta*PI*PI/4.*pow((pow(i/beta,2)+j*j),2)*delta(m,i)*delta(n,j)*(1+delta(0,j))+0.25*gamma*pow((1.+beta*beta),-1.5)*(1.+pow(-1.,j+n))*PI*PI*(pow((i+j),2)*pow((m+n),2)*B1(i+j,m+n)+pow((i+j),2)*pow((m-n),2)*B1(i+j,m-n)+pow((i-j),2)*pow((m+n),2)*B1(i-j,m+n)+pow((i-j),2)*pow((m-n),2)*B1(i-j,m-n));	  
		if (i*j*m*n == 0) {T(p,q) = 0;}
		else {T(p,q) = (1.0*i*n+1.0*i*j*j/n)*B2(m,i)*B2(n,j);}
		
	}
	}

	for (i = 0; i < M; i++)
	{   
	for (j = 0; j < N; j++)
	{
		BC(i*N+j,i) = 1.;
		BC(i*N+j,M+i) = pow(-1,j);
	}
	}

	'''
	extra_code = '''
	float delta(int i, int j)
	{
	if (i == j) {return 1.;}
	else {return 0.;}
	}

	float B1(int i, int j)
	{
	if (i == j) {return 0.5;}
	else {return 0.;}
	}

	float B2(int i, int j)
	{
	if ((i+j)%2==0) {return 0.;}
	else {return i*2./(pow(i,2)-pow(j,2));}
	}

	'''
	weave.inline(code, ['U','T','Nut','M','N','beta','gamma','BC','ds'],
		type_converters=weave.converters.blitz,
		support_code=extra_code,
		headers = ['<math.h>'],
		compiler = 'gcc')

	U = triu(U,1).T+U
	T = triu(T,1).T+T

	K_U = vstack((hstack((U,BC)),hstack((BC.T,zeros((Nbc,Nbc))))))
	K_T = zeros((Nut+Nbc,Nut+Nbc))
	K_T[:Nut,:Nut] = T

	K,A=solve_eigen(K_U,K_T,n_modes)

	return K,A
	
def platebuckling_fixw(M,N,beta,gamma,n_modes=6,ds=0):
	'''
	Elastic Buckling Analysis of four-edge clapped rectangular plate with diagnal stiffeners.
	Input: 
		M,N: Number of series to approximate transverse deflection of the plate;
		beta: aspect ratio of the plate;
		gamma: stiffenning ratio;
		n_modes: Number of buckling modes to calculate;
	Output:
		k: eigen values;
		A: factors of series
	'''

	M,N = int(M),int(N)
	beta,gamma = float(beta),float(gamma)
	Nut = int(M*N)
	U,T = zeros((Nut,Nut)),zeros((Nut,Nut))

	Nbc = 2*M+2*N
	BC = zeros((Nut,Nbc))

	code = '''
	int p,q,i,j,m,n;
	const float PI = 2.*acos(0.); 
	for (p = 0; p < Nut; p++)
	{   
	for (q = p; q < Nut; q++)
	{
		i = p/N; j = p%N;
		m = q/N; n = q%N;
		
		U(p,q) = beta*PI*PI/4.*pow((pow(i/beta,2)+j*j),2)*delta(m,i)*delta(n,j)*(1+delta(i,0)+delta(0,j))+0.25*gamma*pow((1.+beta*beta),-1.5)*(1.+pow(-1.,j+n))*PI*PI*(pow((i+j),2)*pow((m+n),2)*B1(i+j,m+n)+pow((i+j),2)*pow((m-n),2)*B1(i+j,m-n)+pow((i-j),2)*pow((m+n),2)*B1(i-j,m+n)+pow((i-j),2)*pow((m-n),2)*B1(i-j,m-n));	  
		if (i*j*m*n == 0) {T(p,q) = 0;}
		else {T(p,q) = (pow(i,2)*n/m+pow(j,2)*m/n)*B2(m,i)*B2(n,j);}
		
	}
	}

	for (i = 0; i < M; i++)
	{   
	for (j = 0; j < N; j++)
	{
		BC(i*N+j,j) = 1.;
		BC(i*N+j,N+j) = pow(-1,i);
		BC(i*N+j,2*N+i) = 1.;
		BC(i*N+j,2*N+M+i) = pow(-1,j);
	}
	}

	'''
	extra_code = '''
	float delta(int i, int j)
	{
	if (i == j) {return 1.;}
	else {return 0.;}
	}

	float B1(int i, int j)
	{
	if (i == j) {return 0.5;}
	else {return 0.;}
	}

	float B2(int i, int j)
	{
	if ((i+j)%2==0) {return 0.;}
	else {return i*2./(pow(i,2)-pow(j,2));}
	}

	'''
	weave.inline(code, ['U','T','Nut','M','N','beta','gamma','BC','ds'],
		type_converters=weave.converters.blitz,
		support_code=extra_code,
		headers = ['<math.h>'],
		compiler = 'gcc')

	U = triu(U,1).T+U
	T = triu(T,1).T+T

	K_U = vstack((hstack((U,BC)),hstack((BC.T,zeros((Nbc,Nbc))))))
	K_T = zeros((Nut+Nbc,Nut+Nbc))
	K_T[:Nut,:Nut] = T

	K,A=solve_eigen(K_U,K_T,n_modes)

	return K,A
